*
* $Id$
*
* $Log: gntrp.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:39  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:17  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:20:54  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.30  by  S.Giani
*-- Author :
      SUBROUTINE GNTRP (X, PAR, IACT, SNEXT, SNXT, SAFE)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *       COMPUTE DISTANCE UP TO INTERSECTION WITH 'TRAP' VOLUME,  *
C.    *       FROM INSIDE POINT X(1-3) ALONG DIRECTION X(4-6)          *
C.    *                                                                *
C.    *       PAR   (input)  : volume parameters                       *
C.    *       IACT  (input)  : action flag                             *
C.    *         = 0  Compute SAFE only                                 *
C.    *         = 1  Compute SAFE, and SNXT only if SNEXT .GT.new SAFE *
C.    *         = 2  Compute both SAFE and SNXT                        *
C.    *         = 3  Compute SNXT only                                 *
C.    *       SNEXT (input)  : see IACT = 1                            *
C.    *       SNXT  (output) : distance to volume boundary             *
C.    *       SAFE  (output) : shortest distance to any boundary       *
C.    *                                                                *
C.    *    ==>Called by : GNEXT, GTNEXT                                *
C.    *         Author  A.McPherson  *********                         *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gconsp.inc"
      DIMENSION X(6), PAR(11)
#if defined(CERNLIB_SINGLE)
      DIMENSION XT(6), P(11)
#endif
#if !defined(CERNLIB_SINGLE)
      DOUBLE PRECISION XT(6), P(11)
      DOUBLE PRECISION A,B,C,DSN,TOP1,TOP2,BOT1,BOT2
      DOUBLE PRECISION DX1,DX2,X1,X2,X3,DCX1,DCX2,DX,DY
      DOUBLE PRECISION SP1,SP2,DCX,DCY,TANG1,TANG2,A1,B1,A2,B2
      DOUBLE PRECISION XL1,XL2,XL3,XL4,CY,SY,SX,CX,C1,C2
      DOUBLE PRECISION TX,H0,DHDZ,CH
      DOUBLE PRECISION CHTAN,SHTAN
      DOUBLE PRECISION SN1,SN2,SN3,SN4,SN5,SN6,SN7,SN8,P15
#endif
C.
C.               --------------------------------------
C.
      SNXT = BIG
      SAFE = 0.0
      DO 1 I=1,6
   1  XT(I)=X(I)
      DO 2 I=1,11
   2  P(I)=PAR(I)
 
      P15=0.5/P(1)
      DHDZ=(P(4)-P(8))*P15
 
      IF (IACT .LT. 3) THEN
 
 
C       -------------------------------------------------
C       |  Compute safety-distance 'SAFE'  (McPherson)  |
C       -------------------------------------------------
 
 
C            CALCULATE RECTANGLE ON FACE AT Z=-DZ.
 
        DX1=P(4)*P(7)
        X1=-DX1-P(5)
        X2=DX1-P(6)
        IF(X2.GT.X1) X1=X2
        X2=-DX1+P(5)
        X3=DX1+P(6)
        IF(X3.LT.X2) X2=X3
        DCX1=(X1+X2)*0.5
        DX1=(X2-X1)*0.5
C
C            CALCULATE RECTANGLE ON FACE AT Z=DZ.
C
        DX2=P(8)*P(11)
        X1=-DX2-P(9)
        X2=DX2-P(10)
        IF(X2.GT.X1) X1=X2
        X2=-DX2+P(9)
        X3=DX2+P(10)
        IF(X3.LT.X2) X2=X3
        DCX2=(X1+X2)*0.5
        DX2=(X2-X1)*0.5
 
C
        TX=P(2)+(DCX2-DCX1)*P15
C
C            CALCULATE LOCAL RECTANGLE.
C
        SP1=(P(1)-XT(3))*P15
        SP2=(P(1)+XT(3))*P15
C
        DY=P(4)*SP1+P(8)*SP2
        DX=DX1*SP1+DX2*SP2
        DCX=(DCX1+DCX2)*0.5+TX*XT(3)
        DCY=P(3)*XT(3)
C
C            CHECK POINT IS INSIDE RECTANGLE.
C
        IF(ABS(XT(1)-DCX).GT.DX) GO TO 10
C
C            CALCULATE ANGLE OF YZ PLANES.
C
        TANG1=TX+(DX1-DX2)*P15
        TANG2=TX-(DX1-DX2)*P15
C
C            CALCULATE SAFETY FROM YZ PLANES.
C
        SAF1=(XT(1)-DCX+DX)/SQRT(1.0+TANG1*TANG1)
        SAF2=(DX-XT(1)+DCX)/SQRT(1.0+TANG2*TANG2)
C
C            CALCULATE ANGLE OF XZ PLANES.
C
        TANG1=P(3)+DHDZ
        TANG2=P(3)-DHDZ
C
C            CALCULATE SAFETY FROM XZ PLANES.
C
        SAF3=(XT(2)-DCY+DY)/SQRT(1.0+TANG1*TANG1)
        SAF4=(DY-XT(2)+DCY)/SQRT(1.0+TANG2*TANG2)
C
C            CALCULATE SAFETY FROM XY PLANES.
C
        SAF5=P(1)-ABS(XT(3))
C
C            OVERALL SAFETY.
C
        SAFE = MIN(SAF1,SAF2,SAF3,SAF4,SAF5)
 
 
        IF (IACT .EQ. 0) GO TO 999
        IF (IACT .EQ. 1) THEN
          IF (SNEXT .LT. SAFE) GO TO 999
        ENDIF
      ENDIF
 
 
   10 CONTINUE
 
C     ------------------------------------------------
C     |  Compute vector-distance 'SNXT' (McPherson)  |
C     ------------------------------------------------
 
 
 
C           FIRST FIND S TO Z LIMITS.
 
      SN1=BIG
      SN2=BIG
      SN3=BIG
      SN4=BIG
C
      IF(XT(6).EQ.0.0) GOTO 15
      SN1=(P(1)-XT(3))/XT(6)
      SN2=-(P(1)+XT(3))/XT(6)
      IF(SN1.GE.0.) GOTO 15
      SNX=SN1
      SN1=SN2
      SN2=SNX
15    CONTINUE
C
C           NOW Y LIMITS.
C
      H0=(P(4)+P(8))*0.5
C
      TOP1=H0-DHDZ*XT(3)
      TOP2=XT(2)-XT(3)*P(3)
      BOT1=XT(5)-XT(6)*P(3)
      BOT2=DHDZ*XT(6)
C
      IF(BOT1+BOT2.NE.0.0) SN3=(TOP1-TOP2)/(BOT1+BOT2)
      IF(BOT2-BOT1.NE.0.0) SN4=(TOP1+TOP2)/(BOT2-BOT1)
C
C           NOW X LIMTS.
C
      XL1=(P(5)+P(6)+P(9)+P(10))*0.25
      XL2=(P(5)-P(6)+P(9)-P(10))*0.25
      XL3=(P(5)+P(6)-P(9)-P(10))*0.25
      XL4=(P(5)-P(6)-P(9)+P(10))*0.25
C
      CY=XT(2)-XT(3)*P(3)
      SY=XT(5)-XT(6)*P(3)
      CH=H0-DHDZ*XT(3)
C
      A1=XL4*XT(6)*SY+XL3*XT(6)*XT(6)*DHDZ
      B1=XL4*(XT(6)*CY+XT(3)*SY)-XL3*XT(6)*(CH-XT(3)*DHDZ)-XL2*P(1)*SY-
     +XL1*XT(6)*DHDZ*P(1)
      C1=XL4*CY*XT(3)-XL3*XT(3)*CH-XL2*CY*P(1)+XL1*CH*P(1)
C
      CHTAN=(P(7)*P(4)+P(11)*P(8))*0.5
      SHTAN=-(P(7)*P(4)-P(11)*P(8))*P15
      CHTAN=CHTAN+SHTAN*XT(3)
      SHTAN=SHTAN*XT(6)
      CX=XT(1)-XT(3)*P(2)
      SX=XT(4)-XT(6)*P(2)
C
      A2=-P(1)*(DHDZ*XT(6)*SX+SY*SHTAN)
      B2=P(1)*(CH*SX-DHDZ*XT(6)*CX-CY*SHTAN-SY*CHTAN)
      C2=P(1)*(CH*CX-CY*CHTAN)
C
      ISN56=0
      A=A1+A2
      B=B1+B2
      C=C1+C2
      IF(B*B-A*C*4.0.LT.0.0) GO TO 20
      IF(ABS(A).LT.1.0E-7) GO TO 19
      DSN=SQRT(B*B-A*C*4.0)
      SN5=(-B+DSN)*0.5/A
      SN6=(-B-DSN)*0.5/A
      ISN56=1
      GO TO 20
   19 CONTINUE
      IF(ABS(B).LT.1.0E-5) GO TO 20
      ISN56=1
      SN5=-C/B
      SN6=-C/B
   20 CONTINUE
C
      ISN78=0
      A=A1-A2
      B=B1-B2
      C=C1-C2
      IF(B*B-A*C*4.0.LT.0.0) GO TO 30
      DSN=SQRT(B*B-A*C*4.0)
      IF(ABS(A).LT.1.0E-7) GO TO 25
      SN7=(-B+DSN)*0.5/A
      SN8=(-B-DSN)*0.5/A
      ISN78=1
      GO TO 30
   25 CONTINUE
      IF(ABS(B).LT.1.0E-5) GO TO 30
      ISN78=1
      SN7=-C/B
      SN8=-C/B
   30 CONTINUE
C
      IF(SN2.GT.0.0.AND.SN2.LT.SN1) SN1=SN2
      IF(SN3.GT.0.0.AND.SN3.LT.SN1) SN1=SN3
      IF(SN4.GT.0.0.AND.SN4.LT.SN1) SN1=SN4
      IF(ISN56.EQ.0) GO TO 40
      IF(SN5.GT.0.0.AND.SN5.LT.SN1) SN1=SN5
      IF(SN6.GT.0.0.AND.SN6.LT.SN1) SN1=SN6
   40 CONTINUE
      IF(ISN78.EQ.0) GO TO 50
      IF(SN7.GT.0.0.AND.SN7.LT.SN1) SN1=SN7
      IF(SN8.GT.0.0.AND.SN8.LT.SN1) SN1=SN8
   50 CONTINUE
C
      SNXT=SN1
C
  999 CONTINUE
      END
